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Simulations are used to examine the microscopic origins of strain hardening in polymer glasses. 
While traditional entropic network models can be fit to the total stress, their underlying assumptions 
are inconsistent with simulation results. There is a substantial energetic contribution to the stress 
that rises rapidly as segments between entanglements are pulled taut. The thermal component of 
stress is less sensitive to entanglements, mostly irreversible, and directly related to the rate of local 
plastic arrangements. Entangled and unentangled chains show the same strain hardening when 
plotted against the microscopic chain orientation rather than the macroscopic strain. 



The stress needed to deform a polymer glass increases 
as the strain rises. This strain hardening plays a crit- 
ical role in stabilizing polymers against strain localiza- 
tion and fracture, and reduces wear [ij. While models 
have had some success in fitting experimental data, fun- 
damental inconsistencies in fit parameters and trends im- 
ply that our understanding of the microscopic origins of 
strain hardening is far from complete. 

Most theories of strain hardening 0, are based on 
rubber elasticity theory . These entropic network mod- 
els assume that polymer glasses behave like crosslinked 
rubber, with the number of monomers between crosslinks 
equal to the entanglement length N^. The increase in 
the stress a due to deformation by a stretch tensor A 
is attributed to the decrease in entropy as polymers are 

', stretched between affinely displaced entanglements. One 

■ finds 

; a{\)^cjo + GB.9{X)L-\h)/2,h (1) 

where ctq is the yield or plastic flow stress, Gr is the 
' strain hardening modulus, L^^ is the inverse Langevin 
function, g(X) describes the entropy reduction for Gaus- 
sian chains, and L^^{h)/3h corrects for the finite length 
of segments between entanglements. The value of en- 
ters in h, which is the ratio of the Euclidean distance 
between entanglements to the contour length. 

Stress-stretch curves for a wide variety of glassy poly- 
mers can be fit to Eq. [1] but the fit parameters are not 
consistent with the microscopic picture underlying en- 
tropic network models For example, values of 
from fitting h may be several times smaller than those 
, obtained from the plateau modulus G% Q. Entropic 
network models predict Gr w G% near Tg, but mea- 
sured Gr are about 100 times larger Q. Gr also rises 
as T decreases [1, 0] , while any entropic stress must drop 
to zero as T ^ 0. Recent work 0, [1] shows that 
changes in Gr correlate with those in the plastic flow 
stress. Indeed entire stress-stretch curves collapse when 
normalized by (Tq 0] ■ This is not expected from entropic 
models, where ctq is treated as an independent parameter 
arising from local plasticity. A more conceptual difficulty 
in entropic models is that, unlike rubber, glasses are not 
able to dynamically sample chain conformations. 



In this Letter we use simulations to examine the micro- 
scopic origins of strain hardening. While our results for 
the total stress can be fit to Eq. [TJ network models are 
not consistent with observed changes in energy, heat flow 
and molecular conformations. The stress can be divided 
into energetic and thermal contributions. The energetic 
contribution is strictly zero in the entropic model, but 
we find it becomes significant as the segments between 
entanglements are stretched taut. In contrast, entangle- 
ments have little direct influence on the thermal contri- 
bution. This thermal stress is found to be directly related 
to the rate of local plastic rearrangements. Finally, net- 
work models only predict strain hardening for entangled 
chains {N ^ Ne), yet substantial hardening is observed 
for N as small as iVe/4. Results for entangled and unen- 
tangled chains collapse when plotted as a function of the 
microscopic strain-induced orientation of chains rather 
than the macroscopic strain. 

Much of the physics of polymer glasses is independent 
of chemical detail 0, |3, We thus employ a sim- 

ple coarse-grained bead-spring model [ll| that captures 
the key physics of linear homopolymers. Each polymer 
chain contains N beads of mass m. All beads inter- 
act via the truncated and shifted Lennard- Jones poten- 
tial UlAt) = 4uo[(a/r)i2 _ (a/r)'^ - {a/r^f^ + {a/r,)% 
where Vc = 1.5a is the cutoff radius and ULjir) — for 
r > Tc- We express all quantities in terms of the molec- 
ular diameter a, binding energy uq, and characteristic 
time tlj — ^ ma^ /uq. 

Covalent bonds between adjacent monomers on a chain 
are modeled using the finitely extensible nonlinear elastic 
potential U{r) = -{l/2){kRl)ln{l - [r/Rof), with the 
canonical parameter choices Rq = 1.5a and k = 2>Quo/a? 
fll| . Chain stiffness is introduced through a bending po- 
tential Ubend{d) — fcbend(l — cos9), whcrc 9 is the angle 
between consecutive covalent bond vectors along a chain. 
Stiffer chains have lower entanglement lengths. Values 
of Ne obtained from primitive path analysis (PPA) 
range from — 71 for flexible chains {kbend ~ 0) to 

Ne = 22 for kbend = 2.0mo- 

Glassy states were obtained by rapid temperature 
quenches from well-equilibrated melts [l3| to a temper- 
ature T below the glass transition temperature Tg ~ 
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O.SSuoAfcs 1J|- While quench rate affects the initial yield 
stress [l5l |. it had little influence on strain hardening. Pe- 
riodic boundary conditions were imposed, with periods 
Li along directions i = x, y, and z. The initial periods 
were chosen to give zero pressure at T. A Langevin 
thermostat with damping rate 1/tlj was applied to the 
peculiar velocities in all three directions. 

Experiments commonly impose compressive deforma- 
tions because they suppress strain localization 0, [l^, [3] ■ 
Simulations were performed for both uniaxial and plane 
strain compression. Both show the same behavior, and 
only uniaxial results are presented below. The cell is com- 
pressed along z at constant true strain rate e = L^/Lz- 
Results for e — — 10~^/r are shown below, but similar 
behavior is found at e = — 10~^/r. Qualitative changes 
can occur at the higher rates employed in recent atom- 
istic simulations of strain hardening [17 , 1^, 1^, 2^ . The 
stress perpendicular to the compressive axis is main- 
tained at zero by varying and Ly [2l[. Fits to net- 
work models normally assume that the volume remains 
constant and — Ly and this is approximately true in 
our simulations. Then deformation can be expressed in 
terms of a single stretch component A = Lz/L^. 
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FIG. 1: (Color online) Total stress (solid line) and contri- 
butions from heat (dashed line) and potential energy (dot- 
dashed line) for a system with T = 0.2ito/fc_g, A*' — 350 and 
Ne = 22. Fits of a to Eq. [T]with Ne = 13 (dotted line) and 
of CTQ to a straight line are also shown. Both a and g are 
negative under compression. 

Typical strain hardening results are shown in Fig. [T] 
As in experiments, the stress is plotted against .g(A) = 
A^ — 1/A. Then entropic network models (Eq. [T]) attribute 
curvature in the strain hardening regime to reductions 
in entropy associated with the finite length Ne between 
entanglements. The strong upward curvature in Fig. [T] 
can be fit to Eq. [T] (dotted line), but with a value of 
A^e = 13 that is much smaller than that determined from 
G% or PPA (Ne = 22). Similar reductions are found for 
other chain stiffnesses and in fits to experiment. 

The stress represents the incremental work done on a 
unit volume of the system by an incremental strain. It 



can be divided into contributions from changes in the 
internal energy density U and the heat flow out of a unit 
volume Q: a — ajj + aq where ajj — dU/de and ctq — 
dQ jde. The derivation of Eq. [1] assumes that ou does 
not contribute to strain hardening and that oq is entirely 
due to reversible heat associated with changes in entropy. 
Simulations allow these assumptions to be tested. 

Figure [T] shows that results for uq can be fit to the lin- 
ear behavior predicted for the entropy of ideal Gaussian 
chains at \g\ > 1. Fits to smaller \g\ can be obtained with 
A^e = 30 in Eq. [1] but fits to uniaxial and plane strain re- 
sults always give A'e that are larger than values obtained 
from G'^ and PPA, and much larger than values from fits 
to the total stress. Separate simulations show that ctq is 
dominated by irreversible heat flow rather than changes 
in entropy. After straining to a large \g\ the stress is 
returned to zero. The stretch only relaxes about 10% 
and only ~ 5% of the work associated with ctq is recov- 
ered. Similar irreversibility is observed in experiments 
[l6| . confirming that the force can not be entropic. 

The energetic contribution to the stress in Fig. [T] is 
important during the initial rise to the plastic flow stress 
ctq. The value of au then drops to a small constant for 
0.5 < \g\ < 2.5. At larger strains there is a pronounced 
upturn in au that contributes almost all of the curva- 
ture in the total stress. This energetic term thus has a 
crucial effect on fit values of even though the deriva- 
tion of Eq. [1] assumes au = 0. Similar results are found 
for all T and khend, and for uniaxial and plane strain. 
In all cases, aq exhibits nearly ideal Gaussian behavior 
{L~^{h)/3h ~ 1) and au leads to a more rapid rise in 
stress at large stretches. The effect of au increases and 
extends to smaller \g\ as the intrinsic Ne from PPA de- 
creases. 

Examination of the evolving conformations of individ- 
ual chains also provides tests of entropic network mod- 
els. If entanglements act like crosslinks, then polymer 
glasses should deform affinely on scales greater than the 
entanglement spacing. Our recent studies confirm this 
affine displacement, and the associated increase in h as 
segments between entanglements pull taut Fig. [1] 
shows that this increase in h has little effect on the ther- 
mal terms that motivated Eq. [TJ Instead, straightening 
of segments produces large energetic terms by disrupting 
the local packing structure. Energy is stored in increasing 
tension in the covalent bonds countered by compression 
of intermolecular bonds. Experiments also find signifi- 
cant energy storage 3, 22 1, and could in principal track 
au over the full strain hardening regime. 

Further insight into strain hardening is provided by 
examining the dependence on chain length. Entropic 
network models assume that the length should not mat- 
ter for highly entangled systems, N ^ Ne, and there 
should be no network to produce strain hardening for 
N < Ne- Simulations confirm that a is independent of 
A^ for N ^ Ne, but show substantial strain hardening 



3 



for N <NeU, 

Figure [2^ a) illustrates this hardening 
for chains as short as Ne/A. Results for short chains fol- 
low the asymptotic behavior of highly entangled chains 
{N/Nf. > 4) to larger \g\ as N increases. This suggests 
that deformation changes chain conformations on longer 
scales as \g\ increases and that entanglements only be- 
come relevant at large \g\ {\g\ > 2.5 in Fig. [2]). 

The observation of strain hardening implies that the 
microscopic arrangement of chains evolves under stretch- 
ing. One way of quantifying this is through changes 
in the root mean squared components Ri of the end- 
to-end vectors of chains relative to their initial values 
R^. Under an affine deformation at constant volume: 
A = Rz/Rl = {Rl/Rx? = {It/Ru}^. The deforma- 
tion of short chains is subaffine[3, |23|, but we find that 
the above ratios |2J 
five stretch Ag// [25 1 



can all be described by an effec- 
. Figure l^b) shows g{Xeff) as a 
function of g{X) for different N. All chains show sig- 
nificant stretching, and highly entangled systems deform 
affinely. The small deviation between g(A) and g(Ae//) 
for N ^ Ne results from a small increase in density 
4%) at large \g\ rather than nonaffine deformation. 
As N decreases, the deformation becomes subaffine at 
smaller and smaller \g\. This confirms that the scale 
over which chain conformations are distorted grows with 
\g\, and that entanglements only become important for 
\g\ >> 2.5 in this system. 

Strain hardening is directly related to the changes in 
chain conformation represented by Ae/ / rather than the 
macroscopic deformation A [25|. To illustrate this, data 
from Fig. [2jb) are replotted against g{Xeff) in Fig. 
[2jc). Data for different chain lengths collapse onto a 
single curve even though N is as much as 4 times smaller 
than Ne- Similar results are found for other T, and 
for plane strain compression. Deviations are only seen 
when N becomes comparable to the persistence length 
and chains can no longer be viewed as Gaussian random 
walks. These results show that entanglements do not 
have a direct effect on strain hardening. Their main role 
appears to be in forcing the local stretching of chains 
Ae// to follow the global stretch A. 

The recently observed 0, S] correlation between the 
strain hardening modulus Gr and the plastic flow stress 
(To suggests that local plastic rearrangements dissipate 
most of the energy during compression. To monitor the 
rate of plasticity P = Sf/Se, we counted the fraction 
Sf of Lennard- Jones bonds with r < whose length 
changed by more than 20% over small intervals in strain 
5e = 0.005. Tests on this and related amorphous models 
[2^ show that this criterion is large enough to exclude 
elastic deformations, and that 5e is small enough that a 
given bond does not undergo multiple events. To elimi- 
nate plastic rearrangements associated with equilibrium 
aging, the rate of plasticity during deformation was mon- 
itored at T = 0. 

Figure [3] shows the rate of plasticity for Ne — 26 and 





FIG. 2: (Color online) (a) Stress vs. g{\) during uniaxial 
compression at fcsT = .2uo for ktend = 0.75uo, A^e ~ 39 and 
strain rate e — —10'^ /tlj- Successive curves from bottom 

to top are for N = 10 (o), 16 ( ), 25 (A), 40 (■■•), 

70 (squares), 175 (-), and 350 (★). (b) g{\eff) vs. g{X) 
for the same systems. The dot-dashed line corresponds to 
Ae// = A. (c) Stresses for different A'^ collapse when plotted 
against ff(Ae//). 



71. There is a rapid initial rise as a approaches do, fol- 
lowed by a nearly linear increase during the strain hard- 
ening regime. Also plotted in Fig. [3] are results for ctq. 
A fixed vertical rescaling (coincidentally close to unity in 
our units) produces an excellent collapse of P and ctq for 
all Ne- Note that even the fluctuations in the quantities 
are correlated 27]. Similar results arc found for other 
criteria for the rate of plasticity, with only the scaling 
factor changing. 

The above results clearly illustrate that the thermal 
contribution to strain hardening is associated with an 
increase in the rate of plastic rearrangements as chains 
stretch. It remains unclear why this increase should ap- 
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FIG. 3: (Color online) Rate of plastic rearrangements P (solid 
lines) as a function of g{X) for Ne = 26 (upper curve) and 71 
(lower curve) at fcaT — Ouq. Dashed lines show the corre- 
sponding dissipative stress ctq. 



proximately follow the increase in entropic stress pre- 
dicted for Gaussian chains. The entropic stress repre- 
sents the rate of decrease in the logarithm of the num- 
ber of available chain conformations. One possibility is 
that the rate of plastic rearrangements scales in the same 
way because a decrease in the number of conformations 
necessitates larger scale plastic rearrangements. The re- 
lationship between the plastic flow stress and harden- 
ing modulus follows naturally from this picture, and it 
also explains why data for different chain lengths col- 
lapse when plotted against Ae//. Analytic investigations 
of this scenario may prove fruitful. 

Our results for au suggest that the success of Eq. [T] in 
fitting the total stress may be coincidental, and explain 
why fit values of A^e are generally smaller than intrinsic 
values from and PPA. It would be interesting to com- 
pare trends in the calculated ajj and ctq to experimental 
results. These could be obtained by differentiating de- 
formation calorimetry results for the work and heat, but 
existing studie s have only extended to the plastic flow 
regime [Uliilai. 

This material is based upon work supported by the 
National Science Foundation under Grant No. DMR- 
0454947. We thank G. S. Grest for providing equili- 
brated states and K. S. Schweizer and S. S. Sternstein for 
useful discussions. Simulations were performed with the 
LAMMPS software package ( ^http://lammps.sandia.g0Vj ). 
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